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ABSTRACT 

Numerical hydrodynamics calculations are performed to determine conditions under which giant 
planet eccentricities can be excited by parent gas disks. Unlike in other studies, Jupiter-mass planets 
are found to have their eccentricities amplified — provided their orbits start eccentric. We disentangle 
the web of co-rotation, co-orbital, and external resonances to show that this finite-amplitude insta¬ 
bility is consistent with that predicted analytically. Ellipticities can grow until they reach of order 
the disk’s aspect ratio, beyond which the external Lindblad resonances that excite eccentricity are 
weakened by the planet’s increasingly supersonic epicyclic motion. Eorcing the planet to still larger 
eccentricities causes catastrophic eccentricity damping as the planet collides into gap walls. For stan¬ 
dard parameters, the range of eccentricities for instability is modest; the threshold eccentricity for 
growth (^0.04) is not much smaller than the final eccentricity to which orbits grow (^0.07). If this 
threshold eccentricity can be lowered (perhaps by non-barotropic effects), and if the eccentricity driv¬ 
ing documented here survives in 3D, it may robustly explain the low-to-moderate eccentricities 0.1 
exhibited by many giant planets (including Jupiter and Saturn), especially those without planetary 
or stellar companions. 

Subject headings: hydrodynamics — planet-disk interactions — planets and satellites: formation — 
planetary systems: protoplanetary disks 


1. INTRODUCTION 

One of the most surprising revelations of Doppler ex¬ 
oplanet surveys is the prevalence o f Jupiter-mass plan - 
ets on highly elliptical orbits le.g.. iMarcv et akl 120051 1. 
At orbital periods > 10 days, beyond the reach of 
tidal circularization, giant planet eccentricities span 
the full gamut from near-zero to near-unity. There 
is growing evidence that gravitational interactions be¬ 
tween planets can explain the extravagant eccentricities 
observed fe.g.. iT akeda fc Rasiol[20 0^ Juric fc Treirr aind 
[Mil iWu V Lithwickll201~in D awson fc Chiangll2014ll . 

But are planet-planet interactions the whole story? 
After removing observational biases, a substantial frac¬ 
tion of giant planets have low-to-moderate eccentrici¬ 
ties: ~28% have e < 0.05 — our solar system gas gi¬ 
ants belong to this cohort — and fully half have e < 
0.15 (jZakamska et al.l 120111 see their Figure 11, bottom 
panel). T hese statistics are d rawn from the single-planet 
catalog of iButler et al.l l)2006ll . Continued Doppler moni¬ 
toring has not changed the single status of many of these 
planets, particularly at semimajor axes > 1 AU (Bryan 
et ah, submitted). For solitary giants having no stellar 
or planetary perturbers in sight, we look instead to their 
parent gas disks to understand how their ellipticities may 
have arisen. 

Planet-disk interactions are mediated by resonances, 
of which there are as many kinds as there are terms in 
the Fourier expansion of the planet’s potential. Some 
resonances da mp ecc entricity while others excite it. Gol- 
dreich & Sari (120031 hereafter GS03) outlined the cir¬ 
cumstances whereby certain resonant interactions could 
dominate others to excite eccentricity in the net. The 
planet would need to (1) carve out a gap around its or¬ 
bit, and (2) have its eccentricity exceed a threshold value, 
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which GS03 estimated to be on the order of a few per¬ 
cent (see Sections 2 and 4.1 for the technical details). 
Amplification of eccentricity by disk torques could then 
proceed, presumably until the planet crashed into the 
gap walls. Eccentricities excited by disks would then be 
limited by the fractional radial widths of gaps, of order 
^0.1. Planet-disk interactions can thus be argued to be 
relevant for eccentricities in the range ^0.01-0.1. 

The finite-amplitude instability of GS03 has seen little 
if any support in numerical studies. Overwhelmingly, 
planets the mass of Jupiter or lower are seen in nu- 
meric al simulations to have t heir eccentricities damped 
('e.g.. IPanaloizou et al.l 120011: iGresswell fc NelsonI 120061 

Cresswell et al.l 2007: iBitsch fc Klevll20ir)l: iDunhill et al.l 

20131 : iBitsch et al.ll2013fl . Some of these studies found 
that eccentricities grow only for planets of relatively high 
mass > 5-20Mj, via a mechanism that differs from the 
one proposed by GS03 (see Section 2). To our knowl¬ 
edge, t he one numerica l stud y that reported otherwise 
was by iD’Angelo et al.l (|2006fl who found that Jupiter- 
mass planets could have their eccentricities excited to 
values of ~0.1. It is unclear whether their results v indi- 
cate the GS03 mechanism, as ID’Angelo et al.l (|2006ll ob¬ 
served eccentricities to grow starting from zero; in other 
words, no evidence was found for a finite-amplitude in¬ 
stability. 

Many previous numerical studies of eccentricity evolu¬ 
tion used a “live-planet” approach: the planet’s orbit was 
free to evolve under the action of disk torques. Although 
natural enough, a live-planet simulation can be tricky to 
diagnose because all parameters are in flux. We advocate 
here a “fixed-planet” methodology: the planet is kept on 
a fixed eccentric orbit; the disk is allowed to relax to a 
quasi-steady state (one that oscillates consistently with 
the epicyclic phase of the planet); and disk forces on the 
planet are then measured to extract the rate of change 
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of planet eccentricity e as a func tion of e. The fixed - 
pla net approach, as use d by, e.g., iBitsch fc KlevI 1)201011 
and lBitsch et alJ ()2013ll . permits greater control of envi¬ 
ronmental variables and more systematic exploration of 
parameter space. 

In Section[2]we briefly review the theoretical considera¬ 
tions underlying how disks affect planetary eccentricities. 
We summarize the tenets of the theory of GS03 and also 
itemize aspects of the problem that they did not treat. 
Section [3] describes the numerical methods we employed 
to measure e(e) for Jupiter-mass planets. Results, in¬ 
cluding a head-to-head comparison with the predictions 
of GS03, are given in Section S) A summary and outlook 
is contained in Section [5j The Appendix compiles all 
the formulae we used to test GS03, drawn from several 
analytic studies. 

2. THEORETICAL BACKGROUND 

According to GS03, eccentricity excitation requires two 
ingredients: 

• The planet must carve a deep enough gap in the 
disk — or more accurately, a gap with steep enough 
density gradients — that first-order (as expanded 
in the planet’s eccentricity) co-orbital Lindblad 
resonances situated at the very gap center be¬ 
come weaker than first-order external Lindblad res¬ 
onances located roughly a gas scale height h away 
from gap center. The latter drive eccentricity, while 
th e former damp ecce ntricity. Citing calculations 
bv lArtvmowicd l)1993ll . GS03 stated that the gap 
profile must be such that the surface density at the 
locations of the strongest externals must be greater 
than the surface density at gap center by at least 
a factor of ~3 for the externals to defeat the co- 
orbitalsQ 

• First-order co-rotation resonances, which also 
damp eccentricity, must be “saturated” (weak¬ 
ened), meaning that material librating in co¬ 
rotation resonance must not b e replenished by 
visco us inflow of fresh material (jOgilvie fc Lubo-^ 
120031 1. Saturation is effected for sufficiently large 
e > emin; in other words, eccentricity excitation is 
a finite-amplitude instability. 

As estimated by GS03, the minimum eccentricity nec¬ 
essary for e > 0 is 


where E is the disk surface density and v = v/{a?^[j) 
is the kinematic viscosity with dimensions scaled out. 
Equation ([2|) ignores changes in E across the gap which 
can actually be substantial^ we will, in any case, test 
scaling relation ([3|) numerically in Section 4.3. Substi¬ 
tuting dH) into o gives 

emin ~ ■ (4) 


For our standard parameters of g = 0.001 and v = 2.5 x 
10“® (corresponding to a Shakura-Sunyaev a = 0.002 
and disk aspect ratio h/a = 0.036), equation ([T|) — which 
is not meant to be more than an order-of-magnitude es¬ 
timate - gives Cmin "^0.1. 

For e > Cmin, a reasonable expectation not specifically 
discussed by GS03 is that the eccentricity should grow 
until the planet’s radial epicyclic motion causes it to col¬ 
lide with the gap walls. The maximum eccentricity Cmax 
should then scale as w/a0 

A useful order-of-magnitude formula that gives a sense 
of scale is the maximum rate of eccentricit y damping in 
the l imit of small e and no gap clearing ()Artvmowicd 

[ 19 ^ : 

maxle/el^gQ'^flo (5) 


where Eq is the unperturbed disk density. This maxi¬ 
mum rate of eccentricity change is set by the co-orbital 
resonances. Gap clearing can only reduce the magnitude 
of eccentricity changes (and potentially change the sign). 

Apsidal resonances are first-order Lindblad resonances 
with pattern speeds equal to the planet’s apsidal preces¬ 
sion frequency (they hav e wavenumbers m = 1 and £ = 0 
in the Fourier notation of iGoldreich &: Tremainel (11980111. 
They dam p the planet’s eccentricity (|Ward fc HahnI 
119981 l200(Tl l. but are argued by GS03 to be of modest 
importance compared to m > 1 first-order Lindblad res¬ 
onances Q 

Other effects not covered by the linear theory of GS03 
include torques exerted by material in the immediate 
vicinity of the planet, on scales of order the Hill ra¬ 
dius. Circumplanetary material (not necessarily bound 
to the planet) may exert dynamical friction and strongly 
damp the planet’s eccentricity. Properly modeling cir¬ 
cumplanetary flows is challenging and subject to numer¬ 
ical issues such as how the planet’s potential is smoothed 
and how accretion onto the planet is prescribed. An- 



where q is the planet-to-star mass ratio, a is the planet’s 
semimajor axis, Hq is the planet’s orbital frequency, w is 
the gap width, and v is the kinematic viscosity. To es¬ 
timate w, GS03 balance the one-sided principal (zeroth- 
order) Lindblad torque with the local viscous torque 

q^illT,a‘^{a/wf - i^Ea^Ho (2) 

and obtain 

w/a ~ iq^/i>y/^ (3) 

^ We will find in practice that this requirement is met only for 
gaps that are extremely deep in the sense that their central surface 
densities are suppressed by about 3 orders of magnitude relative to 
the background disk. 


^ IFung et al.l (1201^ did account for changes in surface density 
when writing down ((2|l, deriving a scaling relation for gap depth 
that succeeds in reproducing numerical results. These authors re¬ 
placed the left-hand S with Sgap, the right-hand S with the un¬ 
perturbed value Eo, and w with h to arrive at a fairly accurate 
formula for Egap/Eo- Their argument and GSOS’s argument for w 
as presented here are not obviously compatible. 

^ We will present evidence supporting this expectation. Actually, 
we will find that before w/a comes into pla y, th e disk aspect ratio 
h/a < w/a becomes relevant. See Section l4.ll on the supersonic 
weakening of Lindblad resonances and how the weakening leads to 

^max • 

Apsidal (a.k.a. secular) torques will nevertheless be captured 
in our numerical calculations. The disk eccentricity and apsidal 
profile will relax to an equilibrium set by driving from the eccentric 
planet and damping by viscosity; the eccentric disk streamlines will 
backreact secularly onto the planet (probably damping the planet’s 
eccentricity). Although we will not separate out the apsidal/secular 
torque, it is part of t he t o tal t orque that we evaluate from the entire 
disk; see equations ifTsli-ilM. 




































Eccentric Jupiters 


3 


other nonlinear i ssue conc e rns in s tabilities in deep gaps 
(iLi et all l2009t Yu et al.l l2010t iDuffell fc MacFadvenl 


l2013t iFung et al.M20l4 iKanaeawa et alJ 1201511 . If the 

planet mass is large enough or the viscosity is small 
enough, then gap walls can steepen to the point of trig¬ 
gering the Rayleigh instability or the Rossby wave insta¬ 
bility. Gap walls can shed vortices that can stochastically 
torque the planet. 

Finally, we emphasize that the GS03 mechanism for 
eccentricity growth does not align with the common 
view that to drive eccentricity requires near-brown dwarf 
masses and the d ominant influence of the o uter 1:3 Lind- 
blad r esona nce. iPapaloizou et al.l (12001!) . iDunhill et al.l 
(IMl . and iBitsch et al.l (|2013f) found eccentricity driv¬ 
ing only for relatively massive giants (> 5-20Mj); these 
companions opened such wide and deep gaps that they 
interacted primarily with their disks via the outer 1:3 res¬ 
onance, amplifying disk eccentricities which were then 
backreactively sh ared with the planet b y s ecular inter¬ 
action s (see also iKlev fc DirksenI 1200611 . IBitsch et al.l 
(|2013il found eccentricity growth for 5-10 Mj planets 
only when such planets were forced to occupy substan¬ 
tially eccentric orbits, e = 0.2-0.4 (see their Figure 4). 
We will find in the present study that the 1:3 resonance 
is not essential for eccentricity driving; that it is possible 
to excite planetary eccentricities even for Jupiter-mass 
planets, starting with eccentricities as low as a few per¬ 
cent, along the lines envisioned by GS03. 


3. NUMERICAL METHOD 

We address the problem of eccentricity evolution nu¬ 
merically by integrating the 2D (vertically integrated) 
isothermal hydrodynamic equations: 


atS -L V • (EF) = 0 (6) 

dt{'S,Vj) + V • {'Svvj + Pxj — vYiVvj) = —TiVcj) (7) 
P = c^E (8) 


resolution varies with radius to ensure grid cells with 
near-unity aspect ratios, Ar ~ rAcj). 


3.1. Disk Model and Planet Potential 


A simple background disk is employed that ignores gra¬ 
dients in density E, viscosity ly = a{h/aZ, and sound 
speed c: 

II 

o 

W 

II 

(9) 

r2(r) = Do(r/a) 

(10) 

3 

Vr{r) = “2^/^ 

(11) 

P(r) = c^Eo 

(12) 

c = aflo/Ai 

(13) 


where Vr is the background radial accretion velocity and 
A4 = a/h is the constant Mach number, with h the gas 
scale height. 

The gravitational potential at position x is that of the 
star -|- planet: 


cj){x) = GM^ 


^*1 {x — Xp)^ -l- 


(14) 


where q = Mp/M* is the planet-to-star mass ratio and 
e = 0.5h is a smoothing length. The positions of the 
planet Xp{t) and star x,(t) are found by solving Kepler’s 
equation for an eccentric orbit using a Newton-Raphson 
root-finding scheme. Both planet and star are moved 
explicitly in time, keeping the center of mass fixed at 
r = 0. Accretion onto the planet is not modeled. 

Standard model parameters are {q,a,A4} = 

{0.001,0.002,28}. We also vary each of these 3 
parameters separately to values above and below their 
standard value, generating an extra 6 models to explore 
parameter space. 


where E is surface density, P is pressure, v is velocity, 
v is the kinematic viscosity, c is the sound speed, and (j) 
is the gravitational potential from the planet and central 
star. 

The numeri cal integration is carried out us i ng th e 
DISCO code (IDuffell fc MacFadvenl [20Tll \m3. [20ll . 
DISCO is a moving-mesh hydro code that is tailored for 
the study of disks. Computational zones are annular 
wedges that shear past one another to follow the under¬ 
lying flow. By effectively subtracting off the background 
Keplerian flow, DISCO can provide an accurate solution 
for formally supersonic problems, and can integrate for 
long times. 

The numerical domain extends from an inner radius 
Tin = 0.4 to an outer radius Tout = 2, with the planet’s 
semimajor axis located at radius r = a = 10 The domain 
is divided into Nr = 360 logarithmically spaced radial 
zones, corresponding to Ar/r ~ 0.0045. The azimuthal 

® Unless otherwise indicated, we work hereafter in code units: 
GMf = a = 1 (which implies f!o = 1), where the variables have 
their usual meanings. Note that we also set our background surface 
density Eg = 1, which nominally implies a disk mass comparable 
to the stellar mass; but our code ignores self-gravity and therefore 
all of our results for torques, e, and a simply scale as Sq. 


3.2. Caleulating e Numerically 

The planet lives on a fixed eccentric orbit of semimajor 
axis a and eccentricity e (Figure[T]). The code is run until 
the disk surface density relaxes to a pattern that varies 
repeatedly and consistently with the planet’s epicyclic 
motion; typically this takes thousands of planetary or¬ 
bits (see Figure 2). The time derivative of eccentricity is 
time-averaged and recorded as a function of the chosen 
eccentricity, e(e). 

The instantaneous value of e follows from the defini¬ 
tions of the planet’s orbital angular momentum and en¬ 
ergy: 

L = a^DoMpx/l-e^ (15) 

E = -^a'^nlMp. (16) 

Combining the time derivatives of these two quantities 
(and remembering that Dq = \/GM^, /depends on a) 
yields 

e _ P(1 - e2) - 


e 


(17) 
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Fig. 1.— Standard disk-planet system employed in this study: q = 10“®, u = = 2.5x 10“®, A4 = 28 (equivalently, h/a = 0.036). 

The three panels correspond to three choices of planet eccentricity: e = 0.01, 0.05, and 0.12, from left to right. White circles indicate the 
planet’s approximate epicycle. The surface density inside the gap starts as low as Egap/So — 3 x 10“^^ at e = 0.01 and increases with 
increasing e. Eccentricity damps for e = 0.01; amplifies for e = 0.05; and damps for e = 0.12. 
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Fig. 2.— Time derivative of eccentricity as a function of eccentricity for our standard model parameters. Left: e/e as a running average 
over time, demonstrating convergence. Right: Final time-averaged e/e as a function of e. Eccentricity is damped for e < Cmin — 0.04; 
excited for emin < e < emax — 0.07; and damped for e > emax- Thus, there are two attractors: e = 0 and e = emax- Eccentricity damping 
is particularly strong at e > 0.1 when the planet collides into the gap walls. 


where T and P are the torque and power delivered to the 
planet, respectively: 


T = L = rFg 


(18) 


P = E = F-Vp. (19) 

The planet’s velocity is Vp and the disk’s gravitational 
force on the planet is 


F = GM, 




P 

zone j 


(L - 


( 20 ) 


where I is the unit vector pointing from the planet to the 
grid cell of area dAj. The planet’s migration rate can 
also be calculated via 


d _ 2P 
a fl^a^Mp 


( 21 ) 


4. RESULTS 

Residts are first presented for our standard model of 
a Jupiter-mass planet {q = 10“^) in a disk with h/a = 
0.036 {M = 28) and h = v/{a'^nf) = a/M^ = 2.5 x 10"® 
(a = 2 X 10“^). Figure [2] displays the time derivative of 
eccentricity as a function of the eccentricity, e(e). 

The left panel shows the running time average of e/e, 
demonstrating that it can take many thousands of orbits 
to achieve a quasi-steady state (not surprising given the 
low viscosity). The right panel shows the asymptotic 
value of e/e, as a function of e. 

Some highlights from Figure [2j 

• For intermediate eccentricities, e > 0: Jupiter- 
mass planets can, under certain circumstances, 
have their eccentricities excited by the disk. 

• As e —>■ 0, e < 0. Thus e = 0 is an attractor of the 
system for small e. 

• Eccentricity excitation occurs only for e > emin — 
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Fig. 3.— Migration rates as a function of eccentricity for our 
standard model. For modest eccentricities, d is negative and 
roughly independent of e. Once the planet collides with the gap 
walls, d becomes large and positive. The fast outward migration 
coincides with the large negative e seen in Figure [2] 



Eccentricity e 

Fig. 4.— Gap depth Egap/Eo is computed as a function of 
eccentricity. Depths are calculated by averaging E(r) azimuthally 
and over time, excising from the azimuthal average a region of 
radius 0.2a centered on the planet’s guiding center. The minimum 
of this averaged S(r) gives Egap- 

0.04: this is a finite-amplitude instability, as pre¬ 
dicted by GS03. 

• As e increases, eccentricity eventually damps. The 
value Cmax — 0.07 is a second attractor, relevant 
for e > Cinin- 

• For the largest values of e considered, e plunges to 
large negative values. Here the planet’s epicyclic 
motion causes it to collide with the gap walls (see 
also Figure [T]); the gap fills up and eccentricity 
strongly damps. 

Figure [3] plots the migration rate d/a for our stan¬ 
dard model. Note that in contrast with e, the migra¬ 
tion rate does not depend sensitively on e, at least until 
e >w/a 0.1 and the planet crashes into the gap walls. 
The substantial damping of eccentricity for e > w/a 


found in Figure [2] coincides with a large, positive mi¬ 
gration rate in Figure |3 similar to what was observed in 
the live-planet study of lD’ Angelo et al.l 1)200611 . However, 
it should be emphasized that this fast outward migration 
is only sustained as long as the eccentricity is this large. 
In reality the eccentricity should be quickly damped to 
e = Bmax — 0.07 (Figured]), whereupon a < 0 as usual. 

FigureSjshows the gap depth Sgap/So for this system. 
Gap depth is computed by calculating the azimuthally 
averaged and time-averaged surface density as a function 
of radius, S(r), and finding the minimum of this function. 
A region of radius 0.2a centered on the planet’s guiding 
center at (r, 0) = (a, Hgt) is excised from the azimuthal 
average, in order to avoid contamination from material 
very close to the planet. Increasingly eccentric planets 
have shallower gaps. 

In the next section, we elaborate upon all the trends 
highlighted above. We a pply the theory of disk - planet 
interactions pioneered by iGoldreich fc Tremainel (|1980ll 
to see if we can reproduce quantitatiyely the behavior of 
e measured numerically. 

4.1. Detailed Comparison with GS03 
for Standard Model 

Here we compare our numerical results for e for 
our standard model with thos e from analytic the¬ 
ory. U sing formula e derived by IGoldreich k. Tr emainel 
(lIQSOH.lWardl (fTO^.lPapaloizou fc LarwoodI (|2000ll . and 
iQgilvie fc Lubowl (|2003[) . we compute the contributions 
to e from various kinds of resonances: principal Lind- 
blad resonances, first-order (as expanded in the planet’s 
eccentricity) Lindblad resonances, and first-order co¬ 
rotation resonances. Principal co-rotation resonances are 
omitted from our analysis, as these depend on dTi/dr at 
the very gap center; this derivative (difficult to calculate 
reliably) is assumed to be negligibly small for our deep 
gaps. 

The formulae for e are given in the Appendix. They 
depend on surface density I](r) and its slope d'S{r)/dr; 
these two quantities are read directly off snapshots of 
the numerical solution, so in this sense our calculation 
is semi-analytic0 A sampling of surface density profiles 
S(r’) vs. e is provided in Figure [S] overlaid with the lo¬ 
cations of some of the more important resonances. Each 
surface density profile is taken from an individual snap¬ 
shot in time, azimuthally averaged after excising a cir¬ 
cular region of radius = 0.2a centered on the planet’s 
guiding center at {r,9) = (a, Hgt). The excised region 
contains large and highly time-variable overdensities in 
the immediate vicinity of the planet that the analytic 
theories—which govern small disturbances on a smooth 
background—were not intended to treat. We will see 
at the end of this section that torques from this excised 
region are significant in some regions of parameter space. 

The contributions to e from the various resonances are 
dissected in Figure [HI (top panel). As anticipated by 
GS03, the strongest resonances are the first-order Lind¬ 
blad external resonances which excite e, and the first- 
order co-rotation resonances which damp e. The first- 

® In the case of a gapless disk (S = Sq), the equations in the 
.^pendix give a value for e/e < 0 that matches that of equation 
@ to within ~20%, after adjusting the strength of the softening 
term in the generalized Laplace coefficient. 
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Fig. 5. — How gap profiles vary with e for our standard model. 
From bottom to top, surface density profiles correspond to e = 
0.01, 0.03, 0.06, 0.10, and 0.12. Each profile is azimuthally av¬ 
eraged from a late-time snapshot excised of a circular region of 
radius = 0.2a centered on {r,9) = (a, Qq^); the excision removes 
the highly overdense material in the planet’s immediate vicinity 
from the azimuthal average. Resonances from two of the more sig¬ 
nificant wavenumbers (as judged from Figure[7|l are plotted. What 
appear to be nearly overlapping resonances in the figure actually 
do completely overlap in their nominal positions (e.g., the m = 4 
principal Lindblad and m = 4 first-order co-rotation resonances); 
we plot these overlapping resonances with small arbitrary offsets 
for visual clarity only. The co-orbital resonances are so named be¬ 
cause they are located at the planet’s semi-major axis (r = a); they 
should not be confused with the co-rotation resonances, which are 
offset from r = a because they co-rotate with a particular term in 
the planet’s Fourier-expanded potential whose pattern speed does 
not in general equal the planet’s mean motion. 



0.00 0.02 0.04 0.06 0.08 0.10 0.12 
Eccentricity e 


order Lindblad co-orbital resonances also damp e, but 
are weaker because they are situated in the dead center 
of the gap where surface densities are at their lowest. 
Principal Lindblad resonances contribute negligibly to e. 
Figure [7] shows that wavenumbers m ~ 2-8 contribute 
most to e; contributions from higher m, up to our as¬ 
sumed cut-off at TOmax = Af /2, are less important. 

The broad similarity between our semi-analytic calcu¬ 
lation (Figure [51 top panel) and our numerical results 
(Figure ini bottom panel) emboldens us to give the fol¬ 
lowing interpretation of the dynamics. As e increases 
from 0, e switches from negative to positive. This first 
zero crossing is the finite-am plitude instability of GS03 
and lQgilvie fc Lubowl 1)200311 . The instability results be¬ 
cause the first-order co-rotation resonances (which damp 
e) weaken from increasing saturation with increasing e — 
i.e., they weaken following th e F{p ) saturation function 
(see Appendix eauations lA181IA22l plus the discussion at 
the end of this subsection). Above a threshold e, the co¬ 
rotation resonances give way to the first-order Lindblad 
external resonances which render e positive in the net. 
Further increases in e, however, bring e back down to 
a second zero crossing. The external resonances weaken 
as e exceeds h/a, i. e., as the planet’s epicyclic m otion 
becomes supersonic (jPapaloizou fc Larwoodll2000ll . The 
consequence of this supersonic weakening with increasing 
e is that the first-order co-rotation resonances—which 
are unaffected by supersonic motion—together with the 
first-order Lindblad co-orbital resonances regain the up¬ 
per hand at large e to make e < 0. Although the co¬ 
orbital Lindblads suffer from the same supersonic weak¬ 
ening as do the external Lindblads, the co-orbitals yield 


Fig. 6. — Disentangling the web of resonances that contribute 
to eccentricity evolution for our standard disk model. Top: e/e 
as computed semi-analytically from the formulae in the Appendix, 
evaluated using the azimuthally averaged surface density profiles 
S(r) from our numerical calculations (Figure 0 . For e = 0.01, 
we employ Sfr) as computed for e = 0.1, since the surface den¬ 
sity profiles do not change much for e < 0.03. Bottom: e/e com¬ 
puted wholly numerically, with contributions from circumplanetary 
(within 1 Hill radius) and non-circumplanetary material distin¬ 
guished for a few sample e’s. The semi-analytic calculation exhibits 
two trends: (1) a rise in e/e at small e accompanied by a zero cross¬ 
ing that reflects the saturation of first-order co-rotation resonances 
and the growing dominance of first-order Lindblad external reso¬ 
nances; and (2) a drop in e/e at large e accompanied by a second 
zero crossing that reflects the weakening of Lindblad resonances 
from the planet’s increasingly supersonic epicyclic motion. These 
behaviors appear qualitatively reproduced by the numerical calcu¬ 
lations, with modifications introduced by circumplanetary torques 
that linear theory does not capture. The huge drop in e at e > 0.11 
arises from the planet careening into the gap walls. 

a (negative) value of e/e that hardly varies with e; their 
weakening is mitigated by the surface density at gap cen¬ 
ter which grows with e (Figures|3]and|S]), maintaining the 
strength of the co-orbitals. 

Perhaps the most glaring discrepancy between our 
semi-analytic and numerical results is at the largest val¬ 
ues of e. Numerically, at e > 0.11, we find eccentricity 
damping rates that are substantially higher than those 
expected from theory. As the bottom panel of Figure [5] 
indicates, the large negative values of e/e are generated 
from torques exerted by “circumplanetary” material— 
here defined as material within 1 Hill radius of the 
planet’s instantaneous position (the circumplanetary re¬ 
gion so defined is a subset of the excised region used 
to calculate azimuthal averages of surface density). The 
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Fig. 7.— Running sum of e/e vs. azimuthal wavenumber m for 
our standard disk model parameters, calculated semi-analytically. 
Contributions from all kinds of resonances (see Figure[6ll are totaled 
for every m. The sum is truncated at mmax = Af/2 = 14 to 
crude ly account for the “torque cut-ofF (IGoldreich &: Tremaine! 
1198011 . Most of the contributions to e/e arise from m ~ 2—8. The 
locations of the various resonances for m = 4 and m = 7 are shown 
in Figure [5] 

dominance of circumplanetary torques is not surprising 
at large e: the epicyclic motion is so wide that the planet 
collides with the gap walls and suffers dynamical friction 
from high density gas. What is surprising is that circum¬ 
planetary torques also dominate at the smallest value of 
e = 0.01, rendering e more negative than expected from 
theory and pushing the onset of the eccentricity insta¬ 
bility to larger e ~ 0.04 (Figure El bottom panel). The 
properties of the circumplanetary region are uncertain 
and cannot be reliably predicted from linear theory. In 
our numerical calculations, details of the circumplane¬ 
tary flow are subject to such issues as grid resolution, 
smoothing length, and prescriptions for how the planet 
accretes. 

Returning to the first zero-crossing for e as computed 
semi-analytically (top panel of Figure El)! we reiterate 
that it occurs because of co-rotation saturation, as quan¬ 
tified by__F(£)_(egua^n^^^ . This saturation function 
from lOEfilvie fc Lubowl (1200,11 1 (called tc{p) by them) de¬ 
creases with increasing e; it causes the x symbols in the 
top panel of Figure El to approach zero as e increases. 
Physically, F{p) describes how the co-rotation torque 
weakens as viscous diffusion is increasingly unable to 
supply the co-rotation region with fresh librating ma¬ 
terial. The saturation function is separate from the sur¬ 
face density gradient d'E/dr (really vortensity gradient) 
which also facto rs int o the strength of the co-rotation 
torque (equation IA21I) . The surface density gradient at 
the location of a co-rotation resonance also decreases as 
e increases (see Figure El), but the decrease in dE/dr, in 
and of itself, is not as significant as the decrease in F{p). 



0 0.02 0.04 0.06 0.08 0.1 0.12 

Eccentricity e 


Fig. 8.— How e/e varies across parameter space. Data are plot¬ 
ted only for those models whose running time-averages of e/e con¬ 
verged to well-defined values; disks with especially low h/aor low a 
perturbed by planets with high eccentricity exhibited strong insta¬ 
bilities and failed to give convergent answers. Eccentricity driving 
favors large planet masses, small h/a^ and small a, the same region 
of parameter space that produces deep gaps. 

We have shown this by re-computing e/e vs. e using the 
single surface density profile E(r-) evaluated for e = 0.01, 
and obtaining a curve similar to the one shown in the 
top panel of Figure El 

4.2. Dependence on Disk and Planet Parameters 

Figure El shows results for e(e) for our six disk-planet 
systems scattered across parameter space. The observed 
changes to e(e) are complicated and difficult to follow in 
detail. We performed the same semi-analytic analysis for 
thes e models as we did for our standard model (Section 
14.11) , and were able to reproduce only some of the trends 
documented in Figure El Part of our difficulty stemmed 
from circumplanetary torques which often proved signifi¬ 
cant, and which we could not evaluate using the standard 
analytic theory. 

Broadly speaking, however, we can say that e > 0 fa¬ 
vors high-mass planets, thin disks, and low viscosities — 
these cases are highlighted in red in Figure El 








































Duffell & Chiang 


Papaloizou et al (2001) 
d’Angelo et al (2006) 
Bitsch & Kley (2010) 
Dunhillet al (2013) 
Bitsch et al (2013) 
Present Study 


10 ° 10 ^ 10 ^ 10 ° 10 '* 10 ° 10 ° 10 ^ 
K(q,a,h/r) 

Fig. 9. — Several studies of planet-disk eccentricity evolution 
are compared u sing the parameter K = la, which Rov- 

erns gap depths IIDuffell fc MacFadven|[20T3 : IFung et al.ll20l4l1 . A 
rough correlation between high K > 10®-10^ (deep gaps having 
Sgap/So iS 10“®) and e > 0 can be discerned. 

Qualitatively, the circumstances that lead to e > 0 are 
the same ones that produce deep gaps. Gap depths are 
gauged by the dimensionless parameter 

llDuffell &: MacFa.dven 120131 iFung et al.l 12014 
iKanaeawa et al.l 12015 ; iDuffelll I2015D and it is inter¬ 
esting to ask whether this same parameter can predict 
eccentricity growth. Figure [5] lists values of K for 
our seven disk-planet parameter studies, together with 
i^-values from previous studies of eccentricity evolution. 
A rough threshold of AT ^ 10°-10^ dividing e < 0 from 
e > 0 can be discerned — equivalent to a threshold gap 
depth Sgap/So ^ 10“°. 

4.3. Gap Widths 

The width of a gap opened by a planet gives a hard 
upper bound on the eccentricity that can be excited by 
disk torques. As such, it is worthwhile understanding 
how the gap width w depends on input parameters. 

The exercise performed in Section 2 predicts that 
w/a ~ This relation is tested in Figure ITUl 

where gap half-widths are plotted against the dimension¬ 
less parameter /iy. The gap half-width is evaluated by 
differencing the radii at which E = 0.1 Eq and dividing 
by two. The data in Figure [10] appear to conform to a 
power law, but with a somewhat shallower slope than the 
predicted 1/3: w/2a = 0.25((;°/f')°'^^, where w/2 is the 
gap half-width. 

For our standard model parameters, the above fitting 
formula gives w/2a = 0.2. By comparison, the eccentric¬ 
ity beyond which e plummets to large negative values 
(see Figure HI) is e = 0.11; this is a factor of 2 smaller 
than w/2a and suggests that in this context, a more rel¬ 
evant definition for gap half-width would be obtained by 
taking E/Eq = 10“° rather than E/Eq = 10“° — see 
Figure El 

We close with the reminder that the actual value to 
which a planet’s eccentricity relaxes is not given by the 




q^/v 


Fig. 10.— Gap half-widths wf2a (defined as half the distance be¬ 
tween points in the gap where S/Sq = 0.1) increase with increasing 
planet-to-star mass ratio q and decreasing disk viscosity z>. Equa¬ 
tion 0 predicts that ui/2a ~ our measurements are 

fitted by w/2a ~ 0.25{q'^ (solid line). The gap half-width 
gives a hard upper bound on eccentricities that can be sustained 
by planets embedded in their natal gas disks. The actual bound 
Emax on 6 is somewhat smaller — of order a few times h/a — and 
occurs when supersonically-weakened external resonances exactly 
cancel co-orbital and co-rotation resonances to render e = 0. 

gap-collision value, but rather by e„iax (which for our 
standard model equals 0.07) — this is the value for which 
e = 0, and marks where supersonically weakened ex¬ 
ternal Lindblad resonances balance co-rotation and co¬ 
orbital Lindblad resonances fSection Id.lD . 

5. SUMMARY AND DISCUSSION 

This study demonstrates that Jupiter-mass planets 
can have their orbital eccentricities amplified by disk 
torques — provided they open deep enough gaps, and 
provided their eccentricities exceed a threshold value. 
The finite-amplitude instability documented here ap- 
pears to be the same a s that predicted analytically by 
iGoldreich fc Saril (|2003D . Eccentricities are damped by 
first-order co-orbital Lindblad torques and first-order co¬ 
rotation torques. Deep gaps are required to disable 
the former, while finite eccentricities serve to s aturate 
(i.e., weaken) the latter (jOgilvie fc Lubowll2003D . With 
these requirements met, first-order external Lindblad res¬ 
onances can excite a planet’s eccentricitw _ 

Our results are similar to those of iD’Angelo et al.l 
(j200(if ) who also found eccentricity growth for Jovian- 
mass planets at low disk viscosities (a ^ 10“°), but dif¬ 
fer from them insofar as the eccentricity growth that we 
report explicitly requires a non-zero initial eccentricity. 
Other studies did not find eccentricity growth for Jupiter- 
mass planets but used larger viscosities or thicker disks. 
Their results may be reconciled with ours by examining 
the parameter K = q^M^ ja which governs gap depth. 
Eccentricity driving seems to require large K > 10°-10'*, 
i.e., deep gaps of surface density Egap/Eo ^ 1/K. 

Eccentricities do not amplify without bound. As eccen¬ 
tricities increase above the disk aspect ratio h/a — i.e., 
as the planet’s epicyclic motion becomes supersonic — 
the e xternal resonances weaken (iPapaloizou fc LarwoodI 
(200(1) . At the same time, the co-orbital resonances 
strengthen with increasing eccentricity as more disk ma¬ 
terial leaks into the gap. Consequently, eccentricity 
damps above a certain value that scales like h/a] this 







































Eccentric Jupiters 


9 


value is an attractor for the system. In our numerical 
experiments, the attractor eccentricity ranges from 0.07 
to 0.09. At still larger eccentricities — so large that the 
planet collides into gap walls separated by a fractional 
width vjja — the damping of eccentricity becomes catas¬ 
trophically rapid. 

The results of our numerical study align with analytic 
expectations only broadly. Significant torques are ex¬ 
erted by material within a Hill sphere or so of the planet 
that linear theory cannot capture. Modeling circumplan- 
etary flows is technically challenging and we do not claim 
to have gotten it right. In addition to the usual worries 
about smoothing lengths and planetary accretion pre¬ 
scriptions, there looms the possibility that flows in 3D 
could look qualitat ively different from our 2D solutions 
(jOrmel et al.ll2015h . In particular, gaps may be system¬ 
atically s hallower and co-ro tation resonances might never 
saturate (jFung et al.ll2015h . 

The eccentricity driving reported here does not rely 
on non-barotropic (i.e., non-isothermal or non-adiabatic) 
thermodynamics, as our numerical calculations are for 
strictly isothermal disks. Including non-barotropic ef¬ 


fects such as those introduced by external irradiation 
of gap walls may hel p to lower the th reshold eccen¬ 
tricity for instability ([Tsang et al.ll20l3) . If eccentric¬ 
ity driving survives in 3D, it offers a possible expla¬ 
nation for the low-to-moderate eccentricities <0.1 ob¬ 
served for giant planets—including Jupiter and conceiv¬ 
ably Sat urn—without recours e to planet-planet interac¬ 
tions fcf. iTsiganis et al.]l2005D . 
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Planets programs. We are grateful to an anonymous ref¬ 
eree for an encouraging report, and Bertram Bitsch, Re- 
bekah Dawson, Alex Dunhill, Steve Lubow, Alessandro 
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required to compute e. Some of the equations are shared 
between types, but we list the complete set anyway under 
each type for ease of reference. 

First-order Lindblad resonances 
with i = m 1 

The pattern speed of the potential component with 
i = m 1 where m is the a zimuthal wavenumber (see 
[Goldreich fc Tremaine! (|I980f) for their Fourier notation) 
is given by 


n 


m+l,m 


(^) 


n 


planet ■ 


(Al) 


There are two Lindblad resonances associated with this 
potential: a “co-orbital” resonance so called because it 
is located at the planet’s semi-major axis: 


P — C f Uplanet — 11 

^ — ^planet J 


co-orbital 


(A2) 
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and an “external” resonance: 


/3 = 


m — 1 


n = 


m 
m + 1 
TO — 1 


2/3 


- 1 

^planet 


> external. 


(A3) 


We make no accounting for changes in resonance location 
from disk pressure gradients or self-gravity. 

The planet’s eccentricity changes at the rate 


O T'^ 

^ iiplanet planet 


1 


TO (GM,)2 1 + 0.25 [eMf ' 

(A4) 

The second factor invol ving the Mach number A4 = 
Upianet/Zi is taken from iPapaloizou fc LarwoodI 1)2000( 1 
(their equation 32), and accounts for how Lindblad 
torques weaken as the planet’s epicyclic motion becomes 
increasingly supersonic. The first-order Lindblad torque 
equals 




dD 

dr 


,ri 

dr 


-1 


2n 


^771+1,r 




(A5) 


where 


r^^-— = —3fl^ -I- 
dr 


n- 


f^planet , 
TO / 


) ■ (A6) 


First-order Lindblad resonances 
with (. = m — 1 

This case is analogous to the one above. The sum runs 
from TOmin = 2 tO TOmax = Al/2. 


^ ( ~ 1 1 r> 

^^771—1,771 ~ I I ^ ^planet 

m ' 


^ 1 co-orbital 

\L — i ^planet J 


1x2/3 


m 




— I , 1 1 f^planet 

TO -b 1 


> external 


(A9) 

(AlO) 

(All) 


Q rpYi 

^ ^planet planet ?7i —l,?7i ^ 


TO {GM,Y l + 0.25(e7W)^ 

(A12) 


T. 


777— 1,777 


= I r 


d<prn—l^r. 


dD 

dr 


-1 


dr 


dD 

dr 


X 


2^ V 

^ ^777 — 1,777 1 

'*777 — 1,777 / 

(A13) 

TO — 1 ^ 

(AM) 

“^planet I 
m J 


All quantities are evaluated at the resonance location 
(either equation A2 or A3). The surface density E is in¬ 
terpolated from the azimuthally averaged density profile 
of an excised snapshot (see Figure [5] and related text for 
details). The potential amplitude is given by 


/ _ GAfplanet^ 

^777+1,777 — X 


^planet 

- -b TO -b j - 2/3<5mp 


(A7) 


where Sm.i is the Kronecker delta. Note that rdf dr = 
jdd/dp. The Laplace coefficient is 

,m _ 2 r_ cosm(j) _ 

^ ^Jo (1-2/3cos<() + ^2 + 2/7W2)i/2'^‘^ 

(A8) 

and is evaluated numerically. The factor of ac¬ 

counts roughly for how the planet’s poin t-mass potential 
is softened by length ~/i (cf. lWardlll988ll . The coefficient 
of 2 is obtained by calibrating our final answer for e/e 
for the test case of a gapless (E = Eq) disk to match 
approximately the value given by equation ([5]) ; adopting 
a smaller coefficient would overestimate the strength of 
the co-orbital resonances and lead to excessively negative 
values of e. _ 

Equation IA4I is summed from TOmin = 1 to TOmax = 
M. /2, where TOmax crudely appr oximates the torque cut¬ 
off (|Goldreich fc Tremaine! 1980f) . For to = 1 the external 
resonance does not exist. 


^ _ GAfpianet^ 

^777—1,777 — ^ 

^planet 

where is given bv IA81 

First-order Co-rotation Resonances 
with £ = m + I 

The sum runs from TOmin = 1 to TOmax = Ai 12. 
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(A17) 


e = — 


O 

^ ^planet planet '^m+1,777 
TO (GAA)2 


X F(j)) 


(A18) 


The s aturation factor is derived by |Ogilvie_^Ijubow 
(2003) and fitted numerically by IGoldreich fc Sari 
( 1200 ^ : 

(1 + 

(1 + I. 022 p 2)2 ' 
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where 6™2(/3) is given bv IA8I 

First-order Co-rotation Resonances 
with i = m — 1 

The sum runs from = 2 to m^ax = Al/2. 
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Principal Inner Lindblad Resonances 
with £ = m 
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where hy^{l3) is given bv IA8I 

Principal Outer Lindblad Resonances 
with £ = m 

The sum runs from mmin = 2 to mmax = M.I2,. 
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